Linear Regression


In [1]:
from PIL import Image
import numpy as np
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
matplotlib.style.use('bmh')
matplotlib.rcParams['figure.figsize']=(8,5)

簡易的 linear regression 實驗


In [2]:
# 產生隨機數據
X = np.random.normal(0, 3, size=(50,1))
Y = X @ [3] + np.random.normal(0, size=50)
# 畫出來看看
plt.plot(X, Y, 'o');



In [3]:
# 用 numpy 的 lstsq
a = np.linalg.lstsq(X, Y)[0]
a


Out[3]:
array([ 2.95362928])

In [4]:
# 畫出來
plt.plot(X, Y, 'o')
plt.plot(X, X @ a, 'o');


Q

如何加上常數項?

Hint: 使用 np.concatenate, np.ones_like


In [5]:
%run -i q_lstsq.py


a= [  2.96188919  15.04634009]

用 sklearn


In [6]:
from sklearn import linear_model

In [7]:
X = np.random.normal(0, 3, size=(50,1))
Y = X @ [3] + 4 +np.random.normal(0, size=50)

In [8]:
regr = linear_model.LinearRegression()
regr


Out[8]:
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

In [9]:
regr.fit(X,Y)
print(regr.coef_, regr.intercept_)


[ 2.97686497] 4.22689973637

In [10]:
# 畫出來
plt.plot(X, Y, 'o')
plt.plot(X, regr.predict(X), 'o');


Q

畫出 test_X = np.linspace(-10,10, 100) 的圖形


In [11]:
%run -i q_linear_test.py


使用 sklearn 的 datasets


In [12]:
from sklearn import datasets

In [13]:
diabetes = datasets.load_diabetes()
diabetes


Out[13]:
{'data': array([[ 0.03807591,  0.05068012,  0.06169621, ..., -0.00259226,
          0.01990842, -0.01764613],
        [-0.00188202, -0.04464164, -0.05147406, ..., -0.03949338,
         -0.06832974, -0.09220405],
        [ 0.08529891,  0.05068012,  0.04445121, ..., -0.00259226,
          0.00286377, -0.02593034],
        ..., 
        [ 0.04170844,  0.05068012, -0.01590626, ..., -0.01107952,
         -0.04687948,  0.01549073],
        [-0.04547248, -0.04464164,  0.03906215, ...,  0.02655962,
          0.04452837, -0.02593034],
        [-0.04547248, -0.04464164, -0.0730303 , ..., -0.03949338,
         -0.00421986,  0.00306441]]),
 'target': array([ 151.,   75.,  141.,  206.,  135.,   97.,  138.,   63.,  110.,
         310.,  101.,   69.,  179.,  185.,  118.,  171.,  166.,  144.,
          97.,  168.,   68.,   49.,   68.,  245.,  184.,  202.,  137.,
          85.,  131.,  283.,  129.,   59.,  341.,   87.,   65.,  102.,
         265.,  276.,  252.,   90.,  100.,   55.,   61.,   92.,  259.,
          53.,  190.,  142.,   75.,  142.,  155.,  225.,   59.,  104.,
         182.,  128.,   52.,   37.,  170.,  170.,   61.,  144.,   52.,
         128.,   71.,  163.,  150.,   97.,  160.,  178.,   48.,  270.,
         202.,  111.,   85.,   42.,  170.,  200.,  252.,  113.,  143.,
          51.,   52.,  210.,   65.,  141.,   55.,  134.,   42.,  111.,
          98.,  164.,   48.,   96.,   90.,  162.,  150.,  279.,   92.,
          83.,  128.,  102.,  302.,  198.,   95.,   53.,  134.,  144.,
         232.,   81.,  104.,   59.,  246.,  297.,  258.,  229.,  275.,
         281.,  179.,  200.,  200.,  173.,  180.,   84.,  121.,  161.,
          99.,  109.,  115.,  268.,  274.,  158.,  107.,   83.,  103.,
         272.,   85.,  280.,  336.,  281.,  118.,  317.,  235.,   60.,
         174.,  259.,  178.,  128.,   96.,  126.,  288.,   88.,  292.,
          71.,  197.,  186.,   25.,   84.,   96.,  195.,   53.,  217.,
         172.,  131.,  214.,   59.,   70.,  220.,  268.,  152.,   47.,
          74.,  295.,  101.,  151.,  127.,  237.,  225.,   81.,  151.,
         107.,   64.,  138.,  185.,  265.,  101.,  137.,  143.,  141.,
          79.,  292.,  178.,   91.,  116.,   86.,  122.,   72.,  129.,
         142.,   90.,  158.,   39.,  196.,  222.,  277.,   99.,  196.,
         202.,  155.,   77.,  191.,   70.,   73.,   49.,   65.,  263.,
         248.,  296.,  214.,  185.,   78.,   93.,  252.,  150.,   77.,
         208.,   77.,  108.,  160.,   53.,  220.,  154.,  259.,   90.,
         246.,  124.,   67.,   72.,  257.,  262.,  275.,  177.,   71.,
          47.,  187.,  125.,   78.,   51.,  258.,  215.,  303.,  243.,
          91.,  150.,  310.,  153.,  346.,   63.,   89.,   50.,   39.,
         103.,  308.,  116.,  145.,   74.,   45.,  115.,  264.,   87.,
         202.,  127.,  182.,  241.,   66.,   94.,  283.,   64.,  102.,
         200.,  265.,   94.,  230.,  181.,  156.,  233.,   60.,  219.,
          80.,   68.,  332.,  248.,   84.,  200.,   55.,   85.,   89.,
          31.,  129.,   83.,  275.,   65.,  198.,  236.,  253.,  124.,
          44.,  172.,  114.,  142.,  109.,  180.,  144.,  163.,  147.,
          97.,  220.,  190.,  109.,  191.,  122.,  230.,  242.,  248.,
         249.,  192.,  131.,  237.,   78.,  135.,  244.,  199.,  270.,
         164.,   72.,   96.,  306.,   91.,  214.,   95.,  216.,  263.,
         178.,  113.,  200.,  139.,  139.,   88.,  148.,   88.,  243.,
          71.,   77.,  109.,  272.,   60.,   54.,  221.,   90.,  311.,
         281.,  182.,  321.,   58.,  262.,  206.,  233.,  242.,  123.,
         167.,   63.,  197.,   71.,  168.,  140.,  217.,  121.,  235.,
         245.,   40.,   52.,  104.,  132.,   88.,   69.,  219.,   72.,
         201.,  110.,   51.,  277.,   63.,  118.,   69.,  273.,  258.,
          43.,  198.,  242.,  232.,  175.,   93.,  168.,  275.,  293.,
         281.,   72.,  140.,  189.,  181.,  209.,  136.,  261.,  113.,
         131.,  174.,  257.,   55.,   84.,   42.,  146.,  212.,  233.,
          91.,  111.,  152.,  120.,   67.,  310.,   94.,  183.,   66.,
         173.,   72.,   49.,   64.,   48.,  178.,  104.,  132.,  220.,   57.])}

In [14]:
import scipy.stats

In [15]:
scipy.stats.describe(diabetes.target)


Out[15]:
DescribeResult(nobs=442, minmax=(25.0, 346.0), mean=152.13348416289594, variance=5943.331347923785, skewness=0.43906639932477265, kurtosis=-0.8866436055681386)

In [16]:
idx = np.arange(diabetes.data.shape[0])
np.random.shuffle(idx)
X = diabetes.data[idx]
y = diabetes.target[idx]

試試看 linear regression


In [17]:
train_X = X[:-50, 2:3]
train_y = y[:-50]
test_X = X[-50:, 2:3]
test_y = y[-50:]
regr = linear_model.LinearRegression()
regr.fit(train_X, train_y)
plt.plot(train_X, train_y, 'o');
plt.plot(train_X, regr.predict(train_X), 'o');
np.mean((regr.predict(train_X)-train_y)**2)


Out[17]:
4035.1857365613887

In [18]:
plt.plot(test_X, test_y, 'o');
plt.plot(test_X, regr.predict(test_X), 'o');


用所有變數


In [19]:
train_X = X[:-50]
train_y = y[:-50]
test_X = X[-50:]
test_y = y[-50:]
regr = linear_model.LinearRegression()
regr.fit(train_X, train_y)
np.mean((regr.predict(train_X)-train_y)**2)


Out[19]:
2950.6218866898016

In [20]:
np.mean((regr.predict(test_X)-test_y)**2)


Out[20]:
2229.5259631437089

In [21]:
plt.plot(test_X[:, 2:3], test_y, 'o');
plt.plot(test_X[:, 2:3], regr.predict(test_X), 'o');



In [22]:
plt.scatter(regr.predict(train_X), train_y, c='g', s=3)
plt.scatter(regr.predict(test_X), test_y, c='b')
plt.plot([0,300],[0,300],'r', linewidth=1);



In [23]:
groups = np.arange(30,300,60)
predict_y=regr.predict(train_X)
plt.boxplot([train_y[(predict_y>=i-30)&(predict_y< i+30)] for i in groups], labels=groups);
plt.plot(np.arange(1,len(groups)+1), groups,'x');


Regularization

$\frac{1}{2 n} \left\Vert y - Xw\right\Vert_2^2 + α \left\Vert w \right\Vert_1$


In [24]:
regr = linear_model.Lasso(alpha=0.001)
regr.fit(train_X, train_y)
np.mean((regr.predict(train_X)-train_y)**2)


Out[24]:
2950.8234885024399

In [25]:
np.mean((regr.predict(test_X)-test_y)**2)


Out[25]:
2220.756854274488

In [26]:
from sklearn import model_selection
α_space = np.logspace(-4, 0, 50)
scores =[]
for α in α_space:    
    regr.alpha = α
    s = model_selection.cross_val_score(regr, train_X, train_y, cv=3)
    scores.append((s.mean(), s.std()))
scores=np.array(scores).T
plt.semilogx(α_space, scores[0], 'r')
plt.semilogx(α_space, scores[0]+scores[1],'b--')
plt.semilogx(α_space, scores[0]-scores[1],'b--')
plt.fill_between(α_space, scores[0] + scores[1], scores[0] - scores[1], alpha=0.2);



In [27]:
regr = linear_model.LassoCV(alphas = α_space)
regr.fit(train_X, train_y)
print(regr.alpha_)
np.mean((regr.predict(train_X)-train_y)**2)


0.033932217719
Out[27]:
2981.3161667360669

In [28]:
np.mean((regr.predict(test_X)-test_y)**2)


Out[28]:
2126.0973907717112

MNIST

Linear regression 不是個好方法


In [29]:
X = np.random.normal(1, size=(100,1))
y = (X[:,0]>0).ravel()*2-1
regr = linear_model.LinearRegression().fit(X, y)
test_X=np.linspace(-3,3,10).reshape(-1,1)
plt.plot(X, y, 'x');
plt.plot(test_X, regr.predict(test_X), 'r')
plt.plot([-regr.intercept_/regr.coef_[0]]*2, [-1.5,1.5], 'r--')
regr.intercept_


Out[29]:
0.2748505716211514

In [30]:
regr.intercept_


Out[30]:
0.2748505716211514

In [31]:
import gzip
import pickle
with gzip.open('mnist.pkl.gz', 'rb') as f:
    train_set, validation_set, test_set = pickle.load(f, encoding='latin1')
    
train_X, train_y = train_set
test_X, test_y = test_set

In [32]:
regr.fit(train_X, train_y)
regr.predict(test_X)


Out[32]:
array([ 6.27543923,  2.50478365,  1.66828808, ...,  7.80000046,
        5.2956237 ,  5.77673754])

In [33]:
predict_y = np.floor(regr.predict(train_X)+0.5).astype('int').clip(0,9)
np.mean(predict_y == train_y)


Out[33]:
0.23404

In [34]:
predict_y = np.floor(regr.predict(test_X)+0.5).astype('int').clip(0,9)
np.mean(predict_y == test_y)


Out[34]:
0.2286

準確率約 23% 很低

one hot encoding


In [35]:
train_Y = np.zeros(shape=(train_y.shape[0], 10))
train_Y[np.arange(train_y.shape[0]), train_y] = 1

In [36]:
train_y[0]


Out[36]:
5

In [37]:
train_Y[0]


Out[37]:
array([ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.])

In [38]:
from sklearn.preprocessing import OneHotEncoder
onehot_encoder = OneHotEncoder()
onehot_encoder.fit(train_y.reshape(-1,1))
onehot_encoder.transform(train_y.reshape(-1,1)).toarray()[0]


Out[38]:
array([ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.])

In [39]:
# 訓練模型
regr.fit(train_X, train_Y)

# 用 argmax 得到結果
predict_y = np.argmax(regr.predict(train_X), axis=1)
# 計算正確率
np.mean(predict_y == train_y)


Out[39]:
0.85648000000000002

Q

試試看 test accuracy


In [40]:
%run -i q_minst_linear_regression.py


0.8603

Q

用 PCA 試試看